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ABSTRACT 

We performed 3D MHD simulations of planet migration in stratified disks 
using the Godunov code PLUTO, where the disk is turbulent due to the mag- 
netorotational instability. We study the migration for planets with different 
planet-star mass ratios q = Mp/Mg. In agreement with previous studies, for the 
low-mass planet cases (g = 5 x 10^^ and 10~^), migration is dominated by ran- 
dom fluctuations in the torque. For a Jupiter- mass planet (g = Mp/Mg = 10~^ 
for Ms = IMq), we find a reduction of the magnetic stress inside the orbit of 
the planet and around the gap region. After an initial stage where the torque on 
the planet is positive, it reverses and we recover migration rates similar to those 
found in disks where the turbulent viscosity is modelled by an a viscosity. For 
the intermediate-mass planets (g = 5 x 10~^, 10~^ and 2 x 10~^) we find a new 
and so far unexpected behavior. In some cases they experience sustained and 
systematic outwards migration for the entire duration of the simulation. For this 
case, the horseshoe region is resolved and torques coming from the corotation 
region can remain unsaturated due to the stresses in the disk. These stresses are 
generated directly by the magnetic field. The magnitude of the horseshoe drag 
can overcome the negative Lindblad contribution when the local surface density 
profile is flat or increasing outwards, which we see in certain locations in our 
simulations due to the presence of a zonal flow. The intermediate-mass planet is 
migrating radially outwards in locations where there is a positive gradient of a 
pressure bump (zonal flow). 



Subject headings: accretion disks - MHD - planets:formation - turbulence 
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INTRODUCTION 



Understanding why and how fast planets migrate is fundamental to explaining the ob- 
served distribution o f exoplanets and constraining planet formation timescales and efficiencies 
( lAlibert et al.ll2005l ). The basic principle behind the migration of planets in protoplanetary 
disks is the transfer of angular momentum between the planet and its disk. This trans- 
port process occurs at Lind blad resonances and, in a locally isothermal disk, typically leads 
to fa s t inwards migration . ( iGoldreich and Tremaind Il980l : iPapaloizou and Liru Il984j : IWard 
19861 : iTanaka et al.l 120021 ). This is the standard type I migration scenario which applies to 
low- t o intermedia te-mass planets where the specific torque is a linear function of the planet 
mass (]Wardlll997n . 



If the planet is massive enough (the mass depending on the viscosity and the pressure 
scale height), the tidal forces on the disk can eventually overcome the pressure gradient and 
the viscous transport, causing gap opening around the planet orbit. This migration regime, 
referred to as type II, in which the planet-disk interaction can no longer be described by 
linear perturbation theory, is then conceptually very different from the type I regime. Due 
to the gap opening, it is possible for the torques on the planet to cancel in such a way that 
the evolution of the planet's position is determined by the viscous transport o f gas in the 



disk, making the planet move with the disk on viscous timescales (IWardI 119971 : IBate et al. 
2003h . 



Current numerical and analytical calculations estimate migration timescales to be a 
small fraction of the expected disk lifetime, which creates a problem for the survival of plan- 
etary core s. Gas planet core s need to reach a critical mass before the onset of runaway gas 



accretion (jida and Linll2008l ). It is well established that planet population synthesis models 
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survival mechanisms have also been proposed to solve the timescale problem without resort 



ing to an artificially reduce d Type I migration rate (ITerqueml 120031 : iFromang et al.l 12005 
Thommes and Murray! l2006l ) . 



Devia tions from line ar th eory have been found in a numb er of three di r nension al cal- 
culations. iMassetl j^OoJ and b'Angelo et~aD feoosh and later iMasset et all feoOGal ) found 
that for intermediate-mass planets (around Mp = 1 x 10""' M^), the torques on the planet 
can be significantly lower and even reverse sign when the local surface density profile of the 
disk is fiatter (S ~ r" with a = — 0.5) than in the usually assumed Minimum Mass Solar 



Nebula (MMSN) model S 



,-1.5 



. This is found for a certain range of the disk viscosity. In 



this case, the torques from the corotation region can become important. The fiuid elements 
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that are librating (moving in horseshoe orbits in the corotation region) orbit on a U-turn 
trajectory at traihng and leading sides of the planet. These fluid elements exert a torque 
on the planet at each U-turn, which is a symmetric effect on both sides of the planet in an 
inviscid disk; therefore there is no net torque coming from this region after a few librating 
periods. This is referred to as the "saturation" of the corotation torque. In the presence 
of viscosity, if the viscous crossing timescale across the horseshoe region of accreting disk 
material is smaller than the libration timescale, the torques exerted by fluid elements around 
the U-turn are not symmetric at each side of the planet, creating a net positive torque that 
can be sustained. This is refereed to as the "unsaturated" corotation torque, and it can 
depend on the surfa ce density an d on the width of the horseshoe region that in turn depends 
on the planet mass (IWardlll992l ). 

So far most numerical studies of migration and/or gap formation have concentrated on 
the quasi-laminar disk case, where Navier-Stokes shear viscosity is inclu ded in order to mode l 



the viscous stress e s resulting from including turbulence in the disk 
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Nelson et al. (2000); 



Papaloizou and Larwood 



(I2OOOI ): ICrida et al.l ( l2006l ) and many more). There has been strong interest in simulating 
planet-disk interactions in turbulent disks, w here the turbulence is magnetica lly generated 
by the magneto- rotational instability (MRI) (IBalbus and Hawleyl Il99ll . Il998l ). Only ideal 
MHD has been considered so far in global simul ations. The di s k is a ssumed to be fully 
ionized and the magnetic diffusivity is negligible. [Winters et al.l ( 120031 ) looked at gap for- 
mation by intermediate- and large-mass planets in turbulent unstratified disks and the local 
internal stresses around the planet. In the MHD case, they found the gap to be shallower 
and wider compared to the lamina r HP case; the Maxwel l stres ses in the disk dropped in 
the vicinity of the planet's orbit. iPapaloizou and Nelson! (120031 ) performed a comprehen- 
sive study of protoplanets embedded in MHD-turbulent unstratified disks. They found that 
for low mass planets. Type I migration is no longer effective due to large fluctuations in 
the torque. No convergence was reached due to fluctuations of the torque on timescales 
longer than the orbital period and short simulation timescales. However, the torques for 
planets more massive than 30M^ = O.IM 
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Nelson 
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2005) studied the 



long-term evolution of the orbital elements and particularly the excitation of eccentricity 
by turbulent fluctuations. The evolution of the orbital elem ents of particles in MHD turbu- 
lence has also been studied usi ng shearing unstr atified boxes (lYang et al.ll2009l ) and stratified 
boxes including a dead zone ( lOishi et al.l 120071 ). To avoid the expensive MHD simulations, 
other approaches have been taken, such as modeling the tur bulence as a time and space 
varying forcing in a laminar disk model (ILaughlin et al.ll2004[ ). In this case, depending on 
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the amplitude of the forcing, type I migration can be overcome by the random fluctuations 
in the to rque, and randona walk motion will be superimposed on the smooth inward mi- 



gration. iBaruteau and LinI (I2OIOI ) used a similar turbulent forcing model and studied the 



unsaturation of the corotation torque due to turbulence. Depending on the amplitude of the 
turbulence, the corotation torque is found to be unsaturated to a certain level, making the 
total torque increase accordingly (become less negative), slowing down inwards migration. 
Other approaches include the analytical d escription of stochast i c migration of low-mass plan- 
ets using a diffusion-advection equation (jjohnson et al.l l2006l : lAdams and BlochI |2009| ) and 
coupling N-body s imulations with a ra ndom forcin g to study the accre t ion an d formation of 
low- mass planets ( lOgihara et al.l 120071 ). Recently, iNelson and Gressell (120 lOf ) examined the 
velocity dispersion of 1 m to 10 km planetesimals embedded in a turbulent disk, using 3D 
MHD simulations and neglecting stratification, and characterized the stochastic gravitational 
perturbations felt by planetesimals due to MHD turbulence. 

In this paper, we study planet migration in stratified 3D MHD-turbulent disks for planet 
masses in the Type I and II mass range. In Section [2] we describe the numerical setup of 
our simulations, and the initial conditions for the disk and the magnetic field, before the 
addition of a planet. In section 12] we present the results of our simulations and finally in 
section H] we discuss our results. 



2. SIMULATION SETUP 



Simulations whe re performed using the finite volume fluid dynamics code PLUTO 
(IMignone et al.ll2007l ). In the code, time stepping is done using a second order Runge Kutta 
scheme, while the spatial integration is performed using linear interpolation through the 
second order TVD scheme. The Riemann fluxes are computed using the HLLC and HLLD 
solvers for the HD and MHD cases, respectively. The c ode uses the Constrained Transport 
method for preserving a divergence-free magnetic field (IGardiner and Stond 120051) . The nu- 
merical setup for the MHD case follows the setup presented in (IFlock et al.ll2010l ). The MHD 
equations in the isothermal approximation (no energy equation) are given by 



dp „ 



(pv) = 



1 1 
— + (v ■ V)v + -B X (V X B) = — VP - V$„ 

at p p 



v) - (B ■ V)v + (v ■ V)B = 



(1) 
(2) 

(3) 



The potential $g includes contributions from the star and the planet. We work in spher- 
ical coordinates {r,6,(j)), where the computational domain is given by r G [1,10], 6 G 
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[tx/2 - 0.3,7r/2 + 0.3] and G [0,27r]. The grid resolution is {Nr,Ne,N^) = (256, 128,256) 
and it is centered in the center of mass of the planet-star system. The boundary conditions 
for the velocities and magnetic field are periodic in the vertical {6 boundary) and azimuthal 
directions and reflective in the radial direction, except for the transverse magnetic field com- 
ponent, which reverses its sign at the radial boundary. Buffer zones are defined at the radial 
boundaries to avoid boundary effects, where for 1 < r < 2 the magnetic resistivity is given 
by ?7 = 2 X 10"'' (2 — r) and for 9 < r < 10 the resistivity is ?7 = 1 x 10~''(r — 9). 



2.1. Disk and Planet Models 

As an initial condition we take a gas disk in sub-Keplerian rotation around a solar mass 
star. The azimuthal velocity is given by 



v^ = ^vl-c',{a-2b), (4) 

where Vk is the Keplerian velocity and a = 3/2 and b = 0.5 are the exponents of the radial 
power law distribution of the density p oc and sound speed = co(r sin 5)"''. The initial 
density distribution is given by 

p(r, e) = (r sin 9)-^^^ exp (^^^^^^f^^ . (5) 

The disk is described by a locally isothermal equation of state P = c^p. The ratio of the 
pressure scale height h to the radial coordinate of the disk is taken to be a constant such 
that h = H/{r sine) = 0.07. 

The gravitational potential of the planet is given by a softened point-mass potential 

= -Jl 12^ 9M/9 6 

where e is the softening parameter, needed to avoid numerical divergence near the position 
of the planet and 

|rp — = + — 2rpr(sin^pSin^ 

X cos(0p — 0) + cos 9p cos 9) (7) 

is the distance between the planet and a gas particle in the disk. For all the simulations 
e is set to be a fraction of the Hill radius e = krp{Mp/3y^^ with k < 0.5. Table [1] shows 
the parameters of our simulations. Distances are given in units of = lAU, density is 
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given in units of po = 2.6 x lO^^^gcm'^, and velocity is given in units of Keplerian speed 
at lAU, Vq = Vk{lAU). The surface density have been scaled such that the total disk mass 
is O.OlMstar- Magnetic fields are given in units of Bq = ^-np^v^. In the cases where the 
planet is not on a fixed orbit (runs see Tabled]), the equations of motion are integrated with 
a simple leap frog integrator. For the calculation of the torque, we include the entire disk in 
the integration (without Hill sphere tapering, except for simulation R9). The components 
of the torque vector in cartesian coordinates are given by 

= / P^'^ (|r-rpP + e^)3/^ ^^' 

where i G {x, z\ is any of the three cartesian indices and (rp x r)j = (rp x r) ■ Cj, with Cj 
being the cartesian unit vectors. Of course, for studying migration, we are mostly interested 
in the z component of the torque vector. Specific torques are given in units of v\(\AV\ 



2.2. Magnetic Field Configuration 

Before introducing the planet in the simulations, a weak toroidal magnetic field is im- 
posed on the disk given by 

(5„Se,S^) = (0,0,2p/25), (9) 

where p is the initial thermal pressure. This gives an initial azimuthal field with constant 
plasma beta /3 = 25. The field is imposed in a subset of the full computational domain 
given by 2 < r < 9 and 7r/2 — 0.07 < 9 < it /2 + 0.07. The simulation is then followed until 
turbulence generated by the MRI has reached a saturated state. After this stage, we reset 
the density to the initial condition. This is the initial state in which the potential of the 
planet is incorporated and where all our runs start. The azimuthally, vertically and time 
averaged value of the effective a parameter is shown in Figure [1] However, the a stress is 
not constant throughout the vertical dimension. The upper layers of the disk are the most 
active. Figure [2] shows the time evolution of -B^, B'^/Bq and a for simulation RO (see Table 
[1]) that does not include a planet. The top figure shows the characteristic butterfly diagram 
for the azimuthal component of the magnetic field in a turbulent stratified disk. U 



^ A more complete description of the type of model used in this paper and a detailed analysis of the MRI, 
magnetic fields and turbulent spectra can be found in IPlock et al. ( 2011 ). 
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2.2.1. Zonal flows and pressure bumps 



The time-averaged thermal and magnetic pressure and the perturbed (with respect to 
Keplerian) azimuthal velocity are plotted in Figure |3] for run RO. We plot these quantities 
in the mid-plane of the disk, and one scale height above the mid-plane, for a simulation 
without a planet (or, equivalently, a massless planet). The radial gradient of the pressure 
has been removed and the pressure is averaged in the azimuthal direction. As expected 
of zonal flows, we see pressure bumps that correlate wi th bumps in perturbe d azimuthal 
velocity, only phase shifted by one quarter of a period (jjohansen et al.ll2009l ). Bumps in 
thermal pressure correlate with drops in magnetic pressure, a behavior that is seen more 
clearly above the mid-plane, since the MRI is not resolved in the mid-plane of the disk. 
Notice also that in the velocity peaks, the azimuthal velocity exceeds the Keplerian value at 
some radial locations. These structure in the pressure and the velocity lived for the entire 
duration of our simulations, around 1000 inner orbits. These "zonal flows" result from an 
inverse cascade of kinetic energy, e.g. a transport of energy from the MRI unstable medium 
scales, t o the largest scales , whic h is very typica l in a ccretion disks simulations (see for 



instance 



Dzvurkevich et all (120101 1 and iLvra et al.l (l2008f )l 



3. RESULTS 

3.1. Disk torques and migration 

Table [1] summarizes the computational time in local planet orbits for each of the simu- 
lations. The torque was calculated by taking into account the entire disk and its value was 
saved at every time step. We calculated the cumulative average specific torque as 

r„ = ;isLir,At,, (10) 

where r„ is the cumulative average torque up to timestep n and T„ is the total time until 
timestep n. 



3.1.1. Low Mass Planets (q = 5 x 10"'^ and g = 1 x 10 ^) 

Figure H] summarizes the density structure of simulations R2, R5 and R9. Runs Rl 
(g = 5 X 10^®), R2 and R3 {q = 10~^) shows no significant perturbation of the density by 
the planet, and no spiral arms are seen. The turbulent perturbations dominate in this case. 
Figure [5] shows the torque cumulative average torque as a function of local orbits for run 
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Rl. Figure |6] shows the torque for simulations R2 and R3. The fluctuations in the torque 
created by the perturbations in the density can, in both cases, be larger than the mean 
torque expected for standard Type I migration in a laminar disk. Comparing the torque for 
the planet at different positions in the disk, we see that the local (in time) evolution depends 
on the location of the planet. Random variations in the torque can be an order of magnitude 
larger than torques coming from the Lindblad resonances, in addition to the possibility that 
the spiral waves excited at Lindblad resonances are partly or totally suppressed by density 
fluctuations coming from the turbulence, such that the magnitude of this torque can be 
reduced. For the low-mass planet simulations, we find no convergence of the torque on 
timescale of the runs. For a run with a massless planet orbitting at rp = 3.3, a gaussian 
fit of the time distribution of the torque gives a standard deviation of cr ^ 1.5e — 5. We 
also calculate the auto-correlation function of the torque and take the correlation time to 
be given by Tc = /o " ACF{T)dT. This gives Tc ~ 2 local orbits, while the first and second 
zero crossing o f the to r que A CF occur at 0.2 and 0.8 l ocal o rbits. This is in agreement 
with results bvlNelson ( 2005 ) and Fromang and Nelson ( 2009 ) and with estimates used by 
Baruteau and Lin (j2010 ) that are based in previous MHD simulations of turbulent low- mass 



plan et migration. We a l so cal culate the power spectrum of the mid-plane density (see Eq. 



3 in iBaruteau and Liru (120101 )). This is shown in Figure [7] a nd we compare our results 



from MHD simulations to Figure 1 of iBaruteau and LinI ( 120101 ). where the spectrum in the 
result of the forcing model for the turbulence with a ~ 10~^. This comparison is valuable, 
since ultimately, a more complete parameter study migration and turbulence will have to be 
studied in models with forced turbulence. For this a value, in our simulations we find that 
the larger scales carry more power than in the random forcing model, while the two spectrum 
agree for smaller scales for the case that includes the modes with m > 6. The higher power at 
larger scales can result from the higher compressibility of stratified disks, especially at large 
scales, as a vertically stratified disk can respond to compression with vertical expansion. 
However, the overall shape of the spectrum of the MHD simulation agrees better with the 
HD simulation without the m > 6 modes included. Therefore with the proper scaling of the 
amplitude of the turbulence, and a cutoff of these modes, these simulations could reproduce 
the MHD spectrum. Another possibility is that the power at the small scales in the MHD 
simulation is lower due to lacking resolution at these scales. Ultimately, there needs to be a 
physical motivation for the cutoff of the turbulent forcing potential after the first few modes, 
if this is indeed the model that better matches the global MHD simulations. 
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3.1.2. Intermediate Mass Planets (q = 5 x 10 ^, q = 1 x 10"^ andq = 2x IQ-^) 

Figure m shows the log density for run R5. For this simulation (g = 10~^), spiral arms 
are visible and their amplitude is comparable (or larger closer to the planet) to that of the 
perturbations generated by the turbulence. Figure |8] shows the cumulative average specific 
torque for run R4 and Figure [3 shows the torque for simulations R5 and R6. We see that 
in these three simulations there is an initial stage where the torque is negative followed by 
a reversal of the direction of the migration where the torque becomes positive and takes a 
defined value for the rest of the simulation. This happens at different times when we compare 
two different positions of the planet in the disk (runs R5 and R6). Instead of a random walk 
variation in semi major axis superimposed on smooth inwards migration, we find that planets 
of around 30 Earth masses undergo systematic outward migration. This outward migration 
is sustained for the total duration of the simulation. The simulation times for runs R5 and 
R6 are around 600 to 1000 orbits at the inner boundary of the disk (lAU). During this time, 
the density profile in the disk can evolve significantly from the initial state, and although 
the surface density profile can still be fitted by the initial profile (S oc r~^/^), there can be 
changes in the local profile at the position of the planet and accumulation of mass at the 
disk inside the planet's orbit due to turbulent stresses. However, for both simulations, the 
torque is reversed before there is a significant accumulation of mass at the inner boundary 
and it converges to a constant value for the remaining simulation time. The torque for run 
R8 is shown in Figure [TT] In this run the planet mass (g = 2 x 10~^) is now able to modify 
the density profile around its orbit, and opens a partial gap, which affects the convergence of 
the torque. We don't find convergence for the simulation time, but there is still a tendency 
for outwards migration. 

Unlike the simulations for the small-mass planets (Rl, R2 and R3), for the simulations 
R4, R5 and R6, the hill radius of the planet and the horseshoe region are resolved (by ap- 
proximately 4, 4 and 7 grid cells per half width respectively). In this case, the component 
of the torque originating from the horseshoe region can dominate if there is a mechanism 
for keeping the corotation torque unsaturated and the local density profile differs from the 
global profile, possibly increasing outwards, such that the corotation torque can be larger 
than the Lindblad torque, making the total torque positive. There are special locations in 
the disk where is possible for the surface density to increase outwards, due to the appearance 
of zonal flows, as seen in Section 12.2.11 



Com paring the torque valu es of the simulations with analytical estimates by lPaardekooper and Papaloi: 

( 120091 ) or lMasset et al.l ( l2006al ) is not straightforward. First, the undergoing evolution of the 
disk can make the surface density profile at the position of the planet and the effective stress 



- 10 - 



resulting from turbulence vary in time, therefore one torque estimate does not apply at all 
times. On the other hand, the value of the horseshoe drag is very sensitive to the structure 
of the horseshoe region, and the estimate used in the analytical calculations is based on 
a 2D model of the flow around the planet. In our case, the horseshoe region is distorted, 
making the half-width difficult to define. Additionally, the inclusion of ma gnetic fie l ds can 
intr oduce new magnetic re sonances that affect the total torque, as seen by iTerquemI (120031 ) 
and iFromang et al.l (j2005[ ) for a uniform non-turbulent field. For the sake of the compar- 
ison and simplicity, we discard this type of contribution. We attempt a comparison with 
the analytical estimates for the torque, including a contribution of the horseshoe drag. We 
take the total torqu e to be composed of the Lindblad torque in a 3D locally isothermal disk 
dXanaka et aP 120021 ) 



Tund = -(2.340 - 0.099a + 0.4186) (j-^ S^rJU 



4o2 



:ii) 



plus the fully unsaturated non-linear horseshoe dragj ( Paardekooper and Papaloizoull2009l ) 



^Hs = i - <^p<^l- (12) 

Here Qp is the angular frequency of the planet and Sp is the surface density at the 
position of the planet. The cumulative average torque at the end of the simulation for run 
R5 is 2.0 X 10~^ and we take the half width of the horseshoe region to be Xg = 0.25, as is 
measured in our simulations (calculated from the analytical expression, Xg = 0.24). Assuming 
the global surface density profile S oc r~'^, with d = 0.5, will always give a negative torque. 
However, the torque always becomes positive for d = 0.3 and matches the simulation value 
for d = —1.5, which is comparable to the local profile observed in the simulations (see 
middle plot in Figure [T^ . We should also note that already a "close-to-flat" profile can 
significantly reduce the negative torque or change the sign of the torque. For simulation R6, 
the cumulative average torque at the end of the simulation is 3.7 x 10~^ and the measured 
half width of the horseshoe region is Xg = 0.16. In this case, a local profile of d = 0.1 is 
enough to obtain a positive torque, while a local profile of d = —0.4 matches the value of 
the torque obtained in the simulation. However for this simulation, we only observe a fiatter 
profile (see Figure [12]). This discrepancy can be due to the fact that these values are very 
sensitive to the value of Xg, since the horseshoe drag scales as Xg, and we stress that the 
streamlines can be very distorted, therefore making the estimation of Xg difficult. This is 



^We take the expression for an isothermal disk, in the zero gravitational softening limit. 
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a critical parameter, and we find that an increase of 1% to 5% in tlie simulation value of 
Xs with respect to the analytical estimation is enough to reproduce the observed positive 
torques. Therefore, if one assumes that the observed torque is composed of the wave torque 
plus the corotation torque and neglects any additional effect, we see that the corotation 
torque is crucial and able to cancel out or overcome the negative Lindblad contribution for 
standard disk parameters. 

To further test the effect of the local density profile, we performed a simulation with q = 10~^, 
for the planet located at = 4.0 (run R7, Figure ITU]) , initially at the right side of a pressure 
bump (where pressure and density decrease with radius). In this case, the cumulative average 
torque does not clearly converge, and we do not see systematic outwards migration, as the 
cumulative average torque approaches zero. However there is still a significant reduction 
of the torque as compared to the Type I Lindblad torque, which cannot be explained only 
in terms of a locally decreasing radial density profile. This result suggest that even in the 
absence of a pressure bump, inwards migration can be significantly slowed down for this 
planet mass. We note also that we used the expression for the horseshoe drag valid for an 
isothermal disk, so that there is an additional contribution due to the locally isothermal 
profile that we did not take into account. 

To see if the transport of mass in the disk is enough to sustain the unsaturated torque, 
we take the express ion for the minimal a to mantain the unsaturated corotation torque 
dMasset et aPboOGal ) 

am = 0.035g^/2/^-^/^ (13) 

we obtain Om = 0.0003 for q = 10"*^, which is always smaller than what we observe in our 
simulationqj (comparing with the volume average a). For run R6 we also observed that in 
comparison to a purely HD laminar run, in which the planet is able to open a partial gap 
in the disk, the gap in this case is less deep that in the HD case, and also wider, compared 
to the narrower gap seen in the laminar simulation. For run R5, there was no gap opening 
neither in the laminar nor the turbulent runs as the gap opening criterion is not satisfied. We 
also find that the stresses in the disk are affected by the presence of the planet; the volume 
averaged stress decreases as the mass of the planet is increased (see top and middle plots in 
Figure [T3|) . which might be a result of numerical dissipation due to the limited resolution. 



■^However, this expression for am is derived using a 2D model of the HS region, which determines the 
viscous crossing time across the region, the hbration time and the U-turn time. 
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3.1.3. Large Mass Planet (q = 10~^) 

The density structure and spiral arms induced by the Jupiter mass planet in run R9 
{q = 10~^) dominate over the turbulent perturbations and influences the entire disk. Figure 
shows the cumulative average torque. We have excluded the torques coming from the 
Hill sphere for the calculation of the torque, to exclude material bound to the planet which 
is not properly simulated at this resolution or without including other relevant effects. For 
these simulations, initially we see the same trend as for runs R5 and R6, where the torque 
becomes positive, however this is not sustained for this planet mass as torques coming from 
the corotation region are suppressed due to gap opening (see bottom plot in Figure [T^] and 
the right panels in Figure S]). Additionally, the planet modifies the stresses in the disk, and 
therefore, the accretion behavior, as can be seen in the bottom plot in Figure [T3l that shows 
the evolution of the stresses for run R9. In this case, the a stress is progressively suppressed 
and the Reynolds stress dominates the total stress. The reduction of the Maxwell stress is 
seen mostly in the part of the disk inside the planet's orbit and in the gap region (see Figure 
[T^ . Is possible that the magnatic stress is suppressed in the gap region due to the modified 
azimuthal velocity near the planet, which can suppress the MRI locally. Notice that for this 
run, the density has a more laminar appearance (see Figure H]), consistent with a reduction 
of the stress due to the presence of the planet. This could also be a numerical effect that 
appears at this resolution, so further studies at higher resolution are needed. 
In Figure we compare the gap opened by the Jupiter mass planet in a magnetized disk 
(run R9) with an equivalent HD 3D simulation with a viscosity where a = 2 x 10~^ (run 
RIO) and stratification. The time of the snapshots is 100 local orbits. The gap for the 
hydro case is narrower and slightly deeper than the gap formed in the magnetized turbulent 
disk. However the gap is not completely cleaned after this time. We observed the same 
characteristics for low er-mass planets that open only a partial gap in the disk. 



Winters et al.l ( 120031 ) studied a similar case of gap opening, but in a unstratified MHD- 
turbulent disk. In agreement with our results, they found a wider gap when the disk is 
turbulent, and larger transport of mass from the outer to the inner disk (see our Figure [TB]) . 
In contrast to our findings, they find a deeper gap in the hydro case. This can be due to a 
different treatment of the gap opening criteria, since the planet mass in their calculations 
does not satisfy the viscous criterion. In terms of the reduction of the stresses around the 
planet, we find agreement wit h their results. 



Nelson and Papaloizoul ( 120031 ) studied gap opening by a giant planet in an MHD-turbulent 
unstratified disk and compared their results with 2D simulations with an a viscosity. They 
found that a run with an equivalent a stress to the turbulent run produced a shallower 
gap. This is different to our own results where we found the same gap depth. However, we 
should note they only studied how the turbulence affects an already formed gap and did not 
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observe the depletion of the outer disk. The difference to our simulation results could also 
be due to our choice of a in the hydro run. We choose an a matching the global average 
of the turbulent runs, but this is a quantity that varies vertically. In a non-stratified disk 
simulation, this averaging is not necessary. 



4. DISCUSSION AND CONCLUSIONS 

For simulations Rl, R2 and R3, where q = 10~^, during the simulated time, migration 
was dominated by random fluctuations in the torque, that can be orders of magnitude larger 
that what is expected for the value of the Lindb lad or corotat ion torques for this planet 



mass. This is in agreement with simulations by iNelsoru (120051 ) of migration of low-mass 
protoplanets in cylindrical disk models, where stratification is neglected. It is unclear if 
after long term averaging (~ 1000 orbits), the fluctuations will average out to zero while 
some component of the systematic torque will remain. Such a calculation is currently too 
expensive. It will also be difficult to get a steady state without special prescriptions for 
correcting the density, due to the accretion evolution of the disk, in addition to the decrease 
in a stress for long simulation times due to the limited resolution. Another interesting point 
for further studies is to investigate this type of migration with enough resolution to resolve 
the corotation region, to see the impact of the corotation torque in these cases. However, 
even if this torque is present and well resolved, its magnitude would still be small com- 
pared to the amplitude of the fluctuations, since ultimately the torque depends strongly on 
the width of the corotation region, which approaches zero as the planet mass approaches zero. 



As the planet mass is increased by one order of magnitude to g = 10~^, the hill radius 
is now properly resolved and the systematic torque is now large enough to dominate over 
the random component of the torque. Outwards migration in a locally isothermal disk can 
occur due to the viscosity unsaturating the torque coming from the corotation region (where 
the viscous ti mescale across the ho rseshoe region is smaller than the libration timescale), as 
was found by lMasset et al.l (j2006al ) for planets in the intermediate-mass range. Speciflcally, 
a planet with mass ratio q ~ 10"'^ in a disk with h = 0.06 with a flat surface density proflle 
{d = 0) and an a viscosity was found to be the critical mass for which the offset from 
linear theory was the largest. Additionally, a planet mass of q = 10~^ is within the range 
of masses for which reversal of migration occurs, if one extrapolates their results to a disk 
with h = 0.07. However, this is the flrst study that observes the effect of the unsaturated 
corotation torque due to the accretion of mass in the disk provided directly by turbulence that 
has been self-generated by the MRI. For runs R4, R5 and R6, the planet is in locations in the 
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disk where the local surface density profile is either close to flat or increasing outwards, due 
to the pressure bumps seen in Section [2. 2.11 This can make the contribution of the corotation 
torque dominate over the Lindblad torque, which is usually unexpected in an a disk, since for 
realistic density profiles, the Lindblad torque will dominate. This is also consistent with the 
torque in run R6 being initially negative, since at the beginning of the simulation the local 
profile is decreasing outwards, but getting shallower as time increases, eventually reaching 
the point where the torque reverses. For run R5, the slope is almost immediately increasing 
outwards due to the evolution of the disk, which makes the torque positive from the beginning 
of the simulations. We can only roughly compare our numerical results with analytical 
estimates, as was done in the previous section, for the reasons described already there. Also 
in comparing with previous estimates, we also discarded any possible additional contributions 
to the torque that might arise because of the turbulent magnetic fields. The detailed structure 
of the horseshoe region in the presence of turbulence and stratification deserves further study. 
Our results are summarized on Figure [TTl where the torque dependence on planet mass is 
shown. For each simulation, we plot the last value of the cumulative average torque. Note 
however that only for part of the simulations the torque converges to a well defined value. 
It is possible to see a trend of the torque to reverse, corresponding to the addition of the 
contribution of the fully unsaturated horseshoe drag ij^tot = ^Lind + ^hs)- For the plot we 
assumed values for t he width of the horseshoe reg i on th at are 5% larger than the analytical 



estimate given by iPaardekooper and Papaloizod (120091 ) and we use the value of the local 



surface density profile for the calculation of the torque. We see that the trend breaks down 
already for q = 2 x 10~^, where gap opening starts to become important and there is a 
transition into the Type II regime. Error bars represent the standard deviation of the time 
distribution of the torque. Note that since the raw torque is a highly oscillating quantity, the 
standard deviation does not match directly to the amplitude of the turbulent fiuctuations, 
especially in the high-mass planet cases. For run 10, the standard deviation was found to 
be only 20% lower that in the turbulent run R9. In the Type II range, we plot the torque 
corresponding to the viscous timescale of the disk, taking a = 3 x 10"'^. We find reasonable 
agreement with our simulation, taking into account the short simulation time, and that the 
value of the torque is still decreasing in the simulation after 100 orbits. Additionally we use 
the value of the initial, volume averaged a, while the mid-plane value is smaller. 

The question remains about the long term behavior of the torque, and whether this is 
only a transient behavior lasting for the first few hundred orbits (assuming the same local 
surface density profile), afterwards saturating and returning to standard negative Type I 
values. This is still a transient behavior in the sense that the planet can migrate out of the 
part of the disk where the local profile allows for outwards migration and enter a region where 
migration proceeds inwards again. Additionally it is limited by the lifetime of the pressure 
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bumps, which we weren't able to determine. We observe a stable pressure bump through 
the duration of our simulations. If there are other mechanisms such as the ones discussed in 



Masset et al.l (l2006bl ) that produces this type of locally increasing outwards density profile, 



then, in the presence of turbulence, these density bumps can also act as a protoplanet trap 
and halt, slow down or reverse inwards migration. iDzyurkevich et al.l (120101 ) performed non- 
ideal MHD simulations of accretion disks with spatially varying resistivity. They also find 
zonal flows/pressure bumps not only at the snow line, e.g. a region with a jump in resistivity, 
but also inside the more active region. They alr eady suggest tha t smal l planets should get 
trapped at those local pressure maxima (see also iKretke and LinI (120071 )). 



5. SUMMARY 

We studied the migration of planets under the influence of turbulence that is a result 
of the magneto-rotational instability. We find that, under the right conditions, planets can 
undergo systematic outwards migration in a locally isothermal disk. After long term averag- 
ing, transient or long term periods of outwards migration can help the survival and influence 
the mass accretion history of giant planet cores of a certain mass ratio. The contribution 
of the unsaturated horseshoe drag and the stochastic migration of low-mass planets, which 
are both consequences of the turbulence, should be incorporated into planet population syn- 
thesis models in order to test the influence of this element on the produced populations 
of planets. On future work we plan on studying low-mass planet migration in detail using 
similar stratified disk models. 



Giant planets significantly decrease the magnetic stresses in the disk (mostly inside its 
orbit), effectively killing the turbulence, as we observe in our simulations. This is possibly a 
numerical effect and it will affect the accretion behavior of the disk and possibly the Type 
II migration rate of the giant planet. This issue deserves further study, with high resolution 
simulations to determine any possible effects of numerical dissipation of the magnetic fields 
induced by the presence of the planet. Additionally, in agreement with previous studies, we 
find that the gap opened by a planet in the presence of turbulence is wider than the gap 
produced in a quasi-laminar disk with an equivalent a viscosity. 
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Fig. 1. — Initial radial distribution of the time, azimuthally and vertically averaged stress 
parameter (before the addition of the potential of the planet). The dashed and dashed-dot 
lines show the Reynolds Tn^-y and Maxwell Tuax stresses respectively, normalized by the 
initial pressure. The solid line shows the total effective a parameter. 
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Fig. 2. — Time evolution of B^, B^/Bq and a for a run without a planet. The dashed 
and dashed-dotted lines show the Reynolds T^ey and Maxwell TMax stresses respectively, 
normalized to the initial pressure. The solid line shows the total effective a parameter. 
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Fig. 3. — Top figure: Time-averaged thermal and magnetic pressure in the mid-plane (black 
line) and one scale height above the mid-plane (red line). The profiles have been normalized 
to take out the radial variation. Bottom figure: Azimuthal velocity perturbation (with 
respect to Keplerian speed) in the mid-plane (black line) and one scale height above the 
mid-plane (red line). This is from a simulation with no planet included. 



Fig. 4. — Logarithm of the disk density in the mid plane (top row) and in an azimuthal 
cut at the position of the planet (bottom row) for runs R2 (left, q = 10~^), R5 (middle, 
q = 10"^) and R9 (right, q = IQ-^). 




Fig. 5. — Cumulative average torque for run Rl for g = 5 x 10 ^. The red and blue lines 
show the torque exerted by the inner and outer disk respectively. 
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Fig. 6. — Cumulative average torque for runs R2 and R3, for q = 10~^, where the planet is 
located at Vp = 3.3 and = 5.0 respectively. The red and blue lines show the torque exerted 
by the inner and outer disk, respectively. 
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Fig. 7. — Power spectrum of the surface density, averaged in time and azimuthally, from 
MHD simula tion (stars). We co r apare with the power spectrum that results from the turbu- 
lent model of iBaruteau and LinI (120101 ) used in HD simulations, with (triangles) and without 
(crosses) the cutoff of the modes with m > 6, and with effective a ~ 10"^. 
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Fig. 8. — Cumulative average torque for run R4 for g = 5 x 10 ^. The red and blue lines 
show the torque exerted by the inner and outer disk respectively. 
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Fig. 9. — Cumulative average torque for runs R5 and R6, for q = 10~^, where the planet 
is located at Vp = 3.3 and = 5.0 respectively. The red and blue hues show the torque 
exerted by the inner and outer disk, respectively. 
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Fig. 10. — Cumulative average torque for run R7 for q = 10""^, for the planet located at 
Tp = 4.0, initially at the right side of a pressure bump. The red and blue lines show the 
torque exerted by the inner and outer disk respectively. 
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Fig. 11. — Cumulative average torque for run R8 for g = 2 x 10 ^. The red and blue lines 
show the torque exerted by the inner and outer disk respectively. 
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Fig. 12. — Surface density at different times in the simulation. Top, middle and bottom plot 
show the surface density for runs R2, R5 and R9 respectively. The vertical lines shows the 
position of the planet and the extent of the Hill radius. 



-28- 



0.0100 




0.0001 1 , , , , ._ 

20 40 60 80 100 

Planet Orbits 



0.0100 




20 40 60 80 100 

Planet Orbits 



0.0100 E" 




0.0001 1 , , , , ^ 

20 40 60 80 100 

Planet Orbits 



Fig. 13. — Time evolution of the stresses in a disk with an embedded planet. Top, middle and 
bottom plot show the stresses for runs R2, R5 and R9, respectively. The dashed and dashed- 
dotted hues show the Reynolds Tfi^y and Maxwell Tmux stresses, respectively, normahzed to 
the initial pressure. The solid line shows the total effective a parameter. 
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Fig. 14. — Cumulative average torque for run R9 for q = 10^^. The red and blue lines show 
the torque exerted by the inner and outer disk respectively. The torque coming from the 
Hill sphere has been excluded from the calculation. 
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Fig. 15. — Radial distribution of the time, azimuthally and vertically averaged stress pa- 
rameter for run R9. The dashed and dashed-dot lines show the Reynolds T^ey and Maxwell 
Tmux stresses respectively, normalized by the initial pressure. The solid line shows the total 
effective a parameter. 
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Fig. 17. — Specific torque as a function of g = Mp/Mg. Tlie black symbols correspond to 
simulations Rl, R2, R4, R5, R8 and R9, where the position of the planet is Tp = 5.0. The red 
symbols correspond to simulations R3 and R6, where the position of the planet is Vp = 3.3. 
In the g = 5 X 10~^ to g = 2 x 10~^ mass range, we overplot the analytical estimates for the 
torque, taking into account only the Lindblad contribution Ttot = ^Lind (dotted line) and 
both the Lindblad plus the unsaturated horseshoe drag Ttot = ^Lind + ^HS (dash-dotted line), 
for both positions (i.e. local surface density profiles) of the planet, = 5.0 (black line) and 
Tp = 3.3 (red line). For the analytical expressions of the torque, we take the half-width of the 
horseshoe region to be 5% larger than its analytical estimate. The dashed line corresponds 
to the constant Type II migration rate, given by the viscous transport in the disk, using 
a = 2 X 10^'^. Error bars represent the standard deviation of the torque time distribution. 
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Table 1. Simulation Parameters 



Name 


q = Mp/M^ 






Fixed orbit 


Run time (local orbits) 


RO 


Massless 


0.3 


3.3 


Yes 


130 


Rl 


5 X 10-6 


0.3 


5.0 


No 


85 


R2 


10-5 


0.3 


5.0 


No 


140 


R3 


10-5 


0.1 


3.3 


Yes 


89 


R4 


5 X 10-5 


0.3 


5.0 


No 


90 


R5 


10-^ 


0.3 


5.0 


No 


100 


R6 


10-^ 


0.1 


3.3 


Yes 


135 


R7 


10-^ 


0.3 


4.0 


No 


85 


R8 


2 X 10-^ 


0.3 


5.0 


No 


95 


R9 


10-3 


0.3 


5.0 


No 


100 


RIO (HD) 


10-3 


0.3 


5.0 


No 


100 



